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I.  INTRODUCTION 

The  detection  of  nuclear  explosions  in  the  atmosphere 

disturbances^Hencet^tle^research 

extraction  St'iS wlliSeSSS* from  SUrSbS^.phic'sijMls. 

SjSS'»«:  nP?y:"iyd=velop 

digital  dftta  processing  techniques  for  deteFta^  ™  2nd 
analyzing  signals  recorded  by l^^/dTavcf  orms  J{  *collstic. 
^a^r-vSrtr^fafoffi  associated  seismic  surface 
w  ave s  • 

Work  Statement,  with  tasks  a  through  c  designed  to  meet 

°b state-of-the-art6 report^ 'under lining  the1 contributions  of 

research  under  ARPA  to  enhancement  and  understanding  o 
microbarograph  data. 

Prior  to  the  undertaking  of  this  work,  the  recording  o 
microbarograph  data  was  entirely  analog  and  Jr°ceSJg  ®rly 

work  under  the  contract  was  revoiveu  data  without  loss 

of C  informat  ion8  an^digital'programs  were  written  to  dupli- 

as  ^swT^as  22\ts  ^  pjrw&ii 

and  of  direct  digital  processing  techniques. 

The  development  of  diagnostic  parameters  Requires  an 

understanding  of  the  underlying  physics.  Theref or e, 
unaers  tana  mg  studies  were  complemented  with 

waves.  .  ^  nt; 

p  ..,1*,  0f  these  studies  found  application  outsi  e 
this  nroiect  even  Is  the  work  was  going  on  Two  of  the  most 

properties  was  immediately  applicable  to  the  aesig 
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Statement  cf  Work 


Period  of  Performance 


a.  Digitize  microbaragraph  array  data  as 
required  for  research,  and  as  requested 
by  the  Project  Scientist. 

b.  Develop  computer  programs  for  infrasonic 
signal  detection,  analysis,  and  data 
display. 

c.  Evaluate  and  compare  the  various  methods 
for  signal  detection  and  analysis. 

d  Conduct  research  in  the  mechanism  for 
generation  of  Rayleigh  waves  by  nuclear 
explosions  in  the  atmosphere,  with  the 
objective  of  determining  those  character¬ 
istics  of  the  Rayleigh  wave  signal 
spectrum  which  may  be  diagnostic  of  the 
yield  and  height  of  the  explosion. 

e  Conduct  research  in  the  mechanism  for 
generation  of  acoustic/gravity  waves 
in  the  atmosphere  by  nuclear  explosions, 
with  the  objective  of  determining  those 
characteristics  of  the  acous tic/ gravity 
wave  spectrum  which  may  be  diagnostic  ot 
the  yield  r.nd  height  of  the  explosion. 

f.  Survey  and  provide  a  critique  of  contracts 
funded  by  ARPA  in  the  general  field  of 
seismic  and  acoustic  detection  of 
atmospheric  explosions,  as  requested  by 
the  project  officer. 

g.  Conduct  other  research  tasks  as  approved 
by  the  project  officer,  not  to  exceed 
the  estimated  cost  of  the  contract. 


3/1/69  through  2/28/71 


3/1/69  through  2/28/71 


3/1/69  through  2/28/71 
3/1/69  through  3/31/72 


3/1/69  through  3/31/72 


3/1/70  through  3/31/72 


3/1/70  through  2/28/71 
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microbarograph  arrays  being  installed  on  other  projects; 
results  of  f-k  spectral  studies  led  to  a  technique  which, 
when  applied  to  long-period  seismic  data,  permits  us  t° 
strip  off  one  of  two  interfering  signals  to  produce  a  better 
estimate  of  the  second  signal.  Unfortunately,  all  of  ou 
studies  were  not  as  successful.  For  example,  we  failed  in  ou. 
attempts  to  write  a  satisfactory  program  based  „*£*_- 
element  analysis  to  predict  Rayleigh  wave-generation  by  an 
explosion  over  an  atoll.  However,  this  work  is  contmui  g 
a  thesis  program  at  Penn  State  and  should  eventually  be 

successful . 

The  studies  under  this  contract  have  all  resulted  in 
publication,  either  in  technical  journals  or  in  Alexandria 
Laboratory  Reports.  A  list  of  these  publications,  showing 
the  related  tasks  in  the  Work  Statement,  is  given  at  the  end 
of  this  report.  In  the  following  pages,  we  have  outlined 
the  major  results  of  the  program  by  summarizing  the  pertinent 

reports . 
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II.  COMPUTER  PROGRAMS 


Work  on  this  project  included  writing  some  sixty 
computer  programs,  which  were  used  in  our  related  research 
and  are  available  for  use  by  other  research  groups. 

We  have  tabulated  below  brief  summaries  of  the  computer 
programs  developed  by  the  infrasonics  section  of  the 
Alexandria  Laboratories,  or  adopted  from  other  sources 
where  stated.  The  areas  covered  by  these  programs  include 
altering  data  formats,  performing  signal  detection  and 
analysis  in  the  time  domain  and  frequency-wavenumber 
domain,  and  various  theoretical  calculations  involving  the 
acoustic  radiation  from  explosion  sources  in  a  layered 
atmosphere. 

Most  of  the  programs  and  subroutines  are  written  in 
FORTRAN-63  language  for  the  CDC  1604-B  with  some  sub¬ 
routines  coded  in  CODAP,  the  1604  assembly  language.  The 
data  input  is  magnetic  tape  whose  format  is  explained  by 
Kerr  (1971).  The  data  are  output  to  magnetic  tapes.  Output 
displays  are  either  contour  plots  on  the  printer  or  standard 
calcomp  plots. 
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PROGRAMS  AVAILABLE  ON  CONTRACT  NO.  F-14620-69-C-0082 


Info. 

Number 

Program 

Name 

Description 

001 

ESSAFORM 

Reformats  ESSA  array  data  on  VLR  tapes  into  SDL 
unpacked  library  format. 

002 

SLIPSHOD 

Reformats  LASA  slow-mode  infrasonic  data  channels 
into  SDL  format. 

003 

WASBIT 

Digital  simulation  of  ESSA  N-4  correlator  for 
infrasonic  signal  detection  and  analysis. 

004 

FK2D  Tape 

Two  dimensional  wavenumber  spectra  at  fixed 
frequencies . 

005 

VKSPTM2 

Power  as  function  of  frequency  and  wavenumber 
at  fixed  azimuths. 

006 

BREEDER 

Generates  artifical  dnta  consisting  of 
specified  mixtures  of  monochromatic  propa¬ 
gating  plane  waves. 

007 

ARESPONS 

Array  response  in  wavenumber  space. 

008 

ANSWER 

Time- varying  spectrograms,  power  as  a  function 
of  both  frequency  and  time. 

009 

HASH 

Far-field  radiation  from  a  point,  source  in 
earth  atmosphere  system. 

010 

HARKRIDR 

Rayleigh  wave  radiation  from  a  point  pressure 
source  at  the  surface  of  a  layered  earth. 

Oil 

BUMPY 

Two  dimensional  curvature  of  frequency- 
wavenumber  spectra  and  array  response  in 
wavenumber  space. 

012 

CALCULAT 

Predicts  N-4  correlator  output  from  array 
response  in  wavenumber  space. 

013 

PHASDSPN 

Adaptation  of  the  seismic  hyperfine  beam¬ 
steering  program  to  estimate  signal  velocity 
and  azimuth  from  the  correlation  matrix  of 
the  data. 

014 

COHERENCY 

Spectra  matrix  of  multichannel  data  and  display 
of  cross-power  spectra,  coherence  and  phase  as 
functions  of  frequency. 
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Info. 

Number 

Program 

Name 

Description 

015 

PIERCE 

Harkrider's  acoustic-gravity  wave  program  with 
three  dimensional  winds. 

016 

LSQSHIFT 

Hyperfine  beamsteering  parameters  from  data 
correlation  matrix. 

017 

DHCIMERG 

Low-pass  filters,  decimates,  and  merges  two 
consecutive  standard  SDL  subset  tape  records, 
to  facilitate  study  of  long  period  records. 

018 

FKBYBEAM 

Wavenumber  spectra  at  fixed  frequencies  from 
the  beam  sum  variances  of  band-pass  filtered 
data. 

019 

DECODE 

Convert  AFWL  data  into  an  equivalent  1604 

FORTRAN  FORMAT. 

020 

DR-FK2DDS 

Generates  spectral  matrixes  for  any  frequency 
interval,  by  the  Bluestein  algorithm. 

021 

LAPLACE 

Computes  Laplace  transform. 

022 

SAXIER 

Instantaneous  frequency  and  envelope  of  a  time 
series,  by  Cooley-Tukey  algorithm. 

023 

GRANIER 

More  flexible  version  of  Info  004  which  contour 
in  decibels,  centibels  or  milibels  and  permits 
specification  of  wavenumber  region  of  plot. 

024 

LAMB FORM 

Converts  records  from  SDL  format  to  VLR 
format  to  implement  application  of  our  N-4 
correlator  program. 

025 

ANSWER2A 

Plots  power  versus  time  and  frequency  for 
each  of  several  input  channels  and  finally 
plots  the  beam-summed  spectrum. 

026 

DECIMATE 

To  low-pass  filter  and  decimate  a  taped 
digital  record. 

027 

FILTRATE 

To  low-pass,  high-pass,  or  band-pass  filter 
digitally  and  output  a  magnetic  tape  record. 

Four- and  eight-pole  Butterworth  filters. 

028 

DGH83E 

A  faster  variation  of  Info  012  with  facilitated 
input . 

029 

AIRPUNCH 

Calculates  source  height  excitation  functions 
from  eigenfunctions  generated  in  program  Info  010 
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Info. 

Number 

Program 

Name 

Description 

030 

WAYOFF 

Duplicates  airwaves. 

031 

DGHS83 

Generates  seismograms  of  far-field  Rayleigh 
waves  from  surface  pressures  beneath  a  source; 
from  Info  35. 

032 

N4BEAMER 

An  elaboration  of  Info003which  minutely  explores 
the  ESSA  correlator  output  as  it  generates  it, 
and  prints  out  a  bulletin  of  detected  signals, 
with  time,  speed,  bearing,  and  signal- to-noise 
ratio. 

To  generate  an  SDL  subset  magnetic  tape  record 
from  a  multichannel  tape  record  in  alphanumeric 
format  (Isotopes). 

033 

MAGIC 

Calculates  the  Laplacian  of  the  array  response 
at  the  origin  in  wavenumber  space  for  a  given 
array  --  used  in  conjunction  with  Info  013 
(BUMPY) . 

034 

AIRROOTS 

Calculates  eigenvalues  for  normal  modes  of 
acoustic-gravity  waves  in  a  layered  atmosphere. 

035 

MODROOTS 

Calculates  the  real  roots  for  airwaves. 

036 

AIRWAVES 

Calculates  theoretical  pressures  and  velocities 
for  atmospheric  acoustic-gravity  waves  generated 
by  point  sources  in  a  layered  atmosphere. 

037 

PRSHR 

Reads  AFTAC  tapes  of  SHELL  computed  total 
pressures,  reduces  them  to  excess  pressures 
and  writes  them  onto  save  tapes. 

038 

PLOTPRES 

Outputs  the  excess  pressures  from  Info048onto 
the  Calcomp  plotter  as  point  (x)  plots. 

039 

PXFIX 

Reads  the  Inf o 048  save  tapes,  converts  the 
excess  pressures  from  functions  of  variable 
distance  with  time  as  a  parameter  to  functions 
of  variable  time  with  evenly  spaced  distances 
as  parameters  and  writes  them  onto  new  save 
tapes . 

040 

TIMEPRES 

Outputs  the  excess  pressures  from  Info050  onto 
the  Calcomp  plotter  a;;  point  (x)  plots. 
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Info.  Program 

Number  Code  _ Description 


041 

042 

043 

044 

045 

046 

047 

048 

049 

050 

051 


PXTIM 

PLOTPNDT 

PTRANS 


FKFAST 

RESPONSE 

POWER 

MAKETAPE 

DGHS33 

PTRAN52 


Reads  the  Info  050 s ave  tapes,  fits  a  Glasstone 
pulse  to  the  excess  pressure  data,  samples 
the  function  at  regular  time  intervals  and 
writes  them  onto  new  save  tapes. 

Outputs  the  excess  pressures  from  Into  052  onto 
the  Calcomp  plotter  as  line  plots. 

Computes  the  Fourier  transforms  of  the  excess 
pressure  data  on  Info  052  save  tapes  at 
specified  distances  and  frequency  bands  and 
punches  the  resulting  spectra  into  decks  of 
cards . 

A  modification  of  Info  052  for  fitting  a  parabola 
to  the  positive  phase  of  the  pressure  pulse. 

A  modification  of  Info  052  for  fitting  by 
least  squares  a  Glasstone  pulse  to  the  excess 
pressure  data. 

Calculates  frequency-wavenumber  spectra  by  a 
generalization  of  the  process  of  Info  057.  For 
the  case  of  the  13  inner  sensors  at  LAMA, 
this  program  is  more  than  an  order  of  magni¬ 
tude  faster  than  earlier  programs. 

Calculates  a^ray  response  using  the  algorithm 
of  Info058 ,  and  is  also  more  than  an  order  of 
magnitude  faster  for  13  sensors  than  were 
earlier  programs. 

To  compute  power-frequency  spectra;  incorpo¬ 
rates  certain  flexibilities  such  as  control 
of  the  range  of  frequencies  to  plot,  and 
calculates  power  in  units  of  the  time  function 
squared. 

Reformats  specified  data  in  packed  or  unpacked 
library  format,  or  subset  format  into  SDL 
subset  format. 

Modification  of  DGHS83,  to  compute  Rayleigh 
waves  for  a  coupled  nongravitating  ocean- 
solid  earth  half  space. 

Modication  of  PTRANS,  to  compute  the 
pressure  spectra  for  spherically  symmetric 
reflected  Glasstone  pulses. 
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Info.  Program 

Nurobe  r  Code 

052  IBYBHAM 

053  BMSTU 3 

0  54  BURRIlHil; 

0  55  i:m:iBM 

056  IKPOWHR 

0  57  FKSHARCII 

0  58  l-SSASET 

059  BSVI 


Description 


Computes  the  l;-stat  ;  c  as  a  function  oi 
frequency  and  wavenumber  for  multic-.annc 

records . 

Modification  of  existing  program  to  make 
multiple  runs. 

Calculates  the  response  of  a  pipe  array 
by  K.  Burridge,  New  York  Universit,. 

Translates  FORTRAN  languages. 

Spectral  estimation  from  the  best  beams. 
Signal  search  on  three-dimensional  hali- 
spacc  maxima. 

Creates  subset  format  tapes. 

The  program  computes  and  prints  the  vole 
azimuth!  power, "and  Fisher  statists  fo, 

the  best  beam  from  spccilic  inputs. 
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III.  SIGNAL  DETECTION  AND  ANALYSIS 

Out  of  our  evaluation  of  algorithms 

a^o;tantPtools  in  signal  analyses. 

We  describe  four  examples. 

..  n  of  Dhase  velocity  and  apparent  azimuth 
a)  Estimation  ot  pnase  / 

,  r* 1  •  nofiRl  described  a  least- squares 
McCowan  and  ££inn  i  interchannel  time  shifts 

method  for  estimating  the  Ni^  t  piane  wavefront  crossing 
between  corresponding  p  a?®s  technique  made  use 

a  horizontal  array  of  N  elements,  ihe  t  'based  0„  the 

of  estimates  of  the  interchannel  tim< b  s  elements  as  measured 

time  lags  between  all  P°*“5c«relation  matrix.  It  was  shown 
from  the  multichannel  cros' cor  dQm  for  estimation  of 

that  this  led  to  W?,f  «r5Esls  technique  was  then  applied  to 

each  of  the  time  shifts.  Thl is  ‘  1 ’and  apparent  azimuth 

the  estimation  of  the  phase  veio  y  horizontal  array  and 
of  a  plane  wave  propagating: a  ro  McCowan>  1970) :  the 

applied  to  two  5',eI’t*(aJa"kan  earthquake  recoided  on  the 
Ravleigh  waves  from  an  A1,?Mn  e  .J  rray  and  the 
seven-element  UBO  long-period  seismic  ar  7.  in  the 

Icoustic-gravity  «»«nfI"ta0?”£e tiiirtLn-element 
atmosphere,  recorde  P  g0^b  time-domain  and 

microbarographic  array  at  LAMA.  technique  were  used, 

freauercy-domain  formulations  correlation  and  spectral 

i  eTwe  ased  both  the  shifts.  The  results 

matrices  to  measure  the  inter  nrfttion  may  provide  a 

suelett  that  the  time-domain  application  may  p^  velocity 

uslful  algorithm  for  ^nal  detection.  reasPnably  u  with 
dispersion  results  in  both  c as  |  lated  by  another  method, 
theory  and  with  group  velocities  c  eaks  by  linear 

Prewhitening  (flattening  a  means  of  preventing  spectral 
filtering)  was  tested  as  a  “«ans  phase  velocity  and 

leakage  from  ^^^aafs^that  prewhitening  decreases 
azimuth;  the  results  suggest  tPe  ratio  is  high, 

*|»  i^eSSE  when  SS  UgSl-to-noise  ratio  is  low. 

b)  A  high  resolution  f-k  estimator 

We  define  the  curvature  -getrum  ^^derivatives 

Laplacian  of  a  freqUenC^XrndirectiS  is  (McCowan  and  Flinn, 
^Sf.^sSlis  S!.-5S,S:T-!3S  has  the  characteristics 
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of  a  high-resolution  f-k  estimator,  i.e.,  the  spectral  peaks 
are  sharper  than  in  ordinary  f-k  spectra.  The  curvature  of 
an  array  response  function  can  give  additional  information 
about  the  array  capability,  in  addition  to  the  array  response 
function  itself.  In  particular,  defining  the  differential 
array  response  as  the  curvature  of  the  array  respo  isc  at 
k  =  0,  we  show  that  this  differential  response  predicts  the 
angular  resolution  capability  of  an  array.  We  conclude  that 
if  two  arrays  have  differential  responses  I> i  and  1)2,  then 
their  angular  resolution  capabilities  01  and  02  are  related  by 

U101  =  ^*2® 2 ’  Tt  is  shown  that  usin8  the  N'4  correlator  signal 

analysis  algorithm,  the  azimuth  resolution  of  four-  and 
five-element  arrays  with  apertures  as  small  as  20  km  is 
about  one-tenth  of  a  degree.  The  capabilities  of  some  micro- 
barograph  arrays  in  current  use  were  compared. 

c.  Automatic  signal  detection 

A  high  speed  algorithm  for  computation  of  frequency- 
wavenumber  (f-k)  spectra  was  developed  (Smart  and  Flinn, 

1971)  and  two  real-time  infrasonic  data  processing  techniques 
that  it  makes  possible,  were  described:  (1)  Signal  detection 
by  a  search  of  f-k  space.  In  comparison  to  the  N-4  correlator, 
a  broad-band  signal  detector,  the  f-k  search  with  a  Fisher 
detector  has  a  theoretical  advantage,  which  we  verify  in 
practice.  (2)  An  f-k  filter  technique  for  calculating  'best 
beam'  estimates.  This  technique  traces  the  beam  containing 
maximum  power,  from  frequency  to  frequency  through  f-k 
space,  and  thus  allows  for  wandering  of  signal  velocity  and 
arrival  azimuth.  This  maximum  power  function  is  taken  as 
the  frequency  spectrum  cf  the  best  beam.  In  our  programs 
the  Fisher  statistic  of  the  signal  estimate,  and  the  velocity 
and  azimuth,  are  computed  and  displayed  as  functions  of 
frequency.  Examples  from  real  data  for  both  processing 
techniques  were  processed  and  discussed. 

d.  Energy  maxima  in  three  dimensional  f-k  space 

One  of  the  lines  of  research  which  has  had  important 
implications  in  other  geophysical  applications,  concerns 
the  use  of  array  data  to  determine  the  position  of  energy 
maxima  in  three  dimensional  frequency-wave  number  space 
(Smart,  1971).  Algorithms  which  have  evolved  from  this  study 
are  useful  in  separating  two  long  period  signals  which 
traverse  a  seismic  array  at  the  same  time  but  at  only 
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slightly  different  azimuths. 

Research  on  this  topic  was  intensified  when  it  was 
fJ^°Vered  that  two-dimensional  cross-sections  of  finite 
f 1 equency-wavenunber  spectra  can  easily  be  misinterpreted 
since  leakage  of  energy  eccurs  along  lines  of  consM.U 
wavenumber.  In  particular,  the  signal  phase  velocity 
determmed  from  measurements  on  cross-sections  normal  to 
the  frequency  axis  can  be  incorrect. 


The  array  response 
along  lines  of  constant 
function  tends  to  smear 
of  constant  wavenumber, 
velocity  of  a  signal,  a 
first,  a  section  normal 
estimate  the  azimuth  of 


tends  to  smear  the  power  sideways 
period;  the  spectral  smoothing 
the  power  vertically  along  lines 
To  estimate  correctly  the  phase 
two-stage  process  must  be  used: 
to  the  frequency  axis  is  used  to 

section  in  that  a-  th<7  signal;  second,  the  azimuthal 

section  m  that  direction  is  used  to  estimate  signal 

velocity.  In  a  routine  data  processing  procedure8the 

thpTnmnlc section  would  be  formed  at  increments  around 
the  compass.  In  spite  of  the  fact  that  the  combined 
interpretation  of  azimuthal  and  k-plane  sections  suffices 
to  determine  phase  velocity  unambiguously,  this  procedure 

bv^eakaL' !i?WbaC?- :  Wlt£  sPect‘'ai  estimates  contaminated 
y  leakage  along  lines  of  constant  wavenumber  and  planes 

of  constant  period,  this  determination  by  means  of  f-k  and 
k  plane  sections  is  laborious,  We  therefore  devised  a  method 
of  constructing  f reciuency-wavenumber  spectra  in  which  the 

ielocUv  soe?ha?ethCCtrS  alKng  HneS  °{  constant  phase 
in  the  k-pfane  secUon?rU°  PhBSe  Vel°CUy  wU1  be  CVideM 

the  halfwidthaofa?hJ8e  °f  thC  beam'f?rming  algorithm  is  that 
;  naitwidth  of  the  array  response  is  frequency-dependent 

or  eaked  energy.  Thus  if  a  leaked  component  dominates  a 

given  cross-sectien,  its  frequency  will  be  measuraMe  on 

halfwidthSanHCf10"'  The  relaJion  between  apparent  (measured) 
wfdtn  ^  and  frequency  is  that  if  k2  is  the  main  lobe  half 
th  for  a  given  array,  then  the  component  at  f'  will  have 
an  apparent  halfwidth  of  (f0/f»)k2.  Since  f0  frequency 

rl^i-tve-ma)i1KUm  ,a?d  k2  are  known»  the  frequency  f»  may^e 

Width  and  mShlp^^g'by  °f  the  appar°nt  half- 
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IV.  ACOUSTIC  GRAVITY  WAVES 
GENERATED  BY  ATMOSPHERIC  EXPLOSIONS 


a.  Theoretical  models  for  atmospheric  structure 

One  of  the  significant  accomplishments  of  this  project 
was  the  derivation  of  a  new  model  for  atmospheric  structure, 
which  allows  faster  calculation  of  acoustic-gravity  modal 
excitation.  The  standard  programs  (Harkrider,  1964a)  use  a 
model  atmosphere  made  up  of  isothermal  layers,  whereas  the 
new  model  assumes  layers  in  which  the  temperature  varies 
linearly  with  height.  Since  all  atmospheric  structures  contain 
large  intervals  of  height  in  which  the  temperature  varies 
rapidly,  the  isothermal  layer  models  require  many  layers  for 
an  accurate  fit  to  the  structure.  The  modal  excitation  is 
calculated  by  computing  the  product  of  layer  matrices  a  large 
number  of  times,  and  at  each  frequency  several  trigonometric 
functions  must  be  calculated  in  each  layer  matrix.  The 
investment  in  computing  time  is  considerable. 

The  new  method  (Greenfield  and  Harkrider,  1971)  reduces 
the  number  of  layers  required  to  fit  the  structure  by  an 
oder  of  magnitude.  The  layer  matrices  require  calculation 
of  Kummer  functions,  which  require  more  calculation  steps 
than  trigonometric  functions,  but  the  saving  in  tine  and 
storage  space  is  still  impressive. 

In  addition,  the  new  formulation  allows  precise  study 
of  Tolstoy's  (1967)  prediction  that  strong  positive  tempera¬ 
ture  gradients  should  cause  anomalous  acoustic-gravity  wave 
propagation,  in  that  the  acoustic  cutoff  period  (the  short- 
period  cutoff)  is  greater  than  the  Brunt  Vaisala  period. 

This  study  is  being  continued  at  the  California  Institute 
of  Technology  under  another  contract. 

b.  Atmospheric  waves  as  functions  of  yield  and  height  of  burst 

The  continuing  testing  of  nuclear  weapons  in  the  atmos¬ 
phere  has  stimulated  interest  in  the  use  of  infrasonic  data 
for  diagnosing  yields  and  burst  altitudes.  We  calculated 
theoretical  barograms  for  a  variety  of  yields  and  burst 
heights,  using  programs  supplied  by  D.G.  Harkrider  (1964a). 

Harkrider's  programs  were  modified  to  include  the  energy 
injection  source  described  by  Pierce  (1968).  The  programs  are 
listed  in  an  appendix  of  the  detailed  report  by  Newton 
et  al.  (1972). 
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Using  the  energy  injection  source,  theoretical  spectra 
were  computed  for  point  sources  in  an  ARDC  model  atmosphere 
(Wares  et  al.,  1960),  for  an  epicentral  distance  of  10,000  km. 
Excitation  was  computed  for  the  modes  which  propagate  in  the 
period  range  of  the  microbarograph  stations,  0.4  to  14  minutes 
these  include  the  gravity  mode  GRq  and  acoustic  modes  So  to 
Sg.  The  individual  modes  were  summed  to  yield  composite 
barograms . 

We  found  that  th  long-period  power  (greater  than 
1  minute  period)  varies  as  the  square  of  the  yield  for  low- 
altitude  explosions,  and  as  the  two-thirds  power  of  yield 
for  high  altitudes  (around  100  km).  For  a  given  yield,  the 
long-period  power  increases  with  increasing  burst  height  to 
a  maximum  and  then  decreases;  the  maximum  occurs  at  lower 
altitudes  for  higher  yields  (Newton  et  al.,  1972). 


c.  Spacial  coherency  of  acoustic  gravity  waves 

The  coherence  of  atmospheric  acoustic-gravity  waves 
has  been  measured  in  the  period  range  10-100s  at  the  Large 
Aperture  Microbarograph  Array  in  south-eastern  Montana.  The 
acoustic- gravity  waves  observed  were  signals  generated  by 
presumed  nuclear  explosions.  The  decrease  of  coherence 
with  increasing  distance  between  pairs  of  microbarographs 
is  less  rapid  in  the  direction  of  wave  propagation  than 
transverse  to  it.  Variation  of  direction  of  arriyal  over 
a  small  range  of  azimuth  (+5°)  explains  the  spatial 
behaviour  of  coherence  in  the  direction  normal  to  the 
wave  propagation;  variation  of  phase  velocity  of  +10m/sec 
explains  the  behaviour  along  the  direction  of  wave  propaga¬ 
tion.  Both  effects  may  be  due  to  inhomogeneities  in  the 
atmosphere;  the  velocity  variation  may  be  due  to  the  presence 
in  the  signal  of  several  normal  modes  of  acoustic-gravity 
waves,  each  travelling  at  a  slightly  different  phase  velocity 
in  the  range  300-330  m/sec. 
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V.  RAYLEIGH  WAVES  GENFRATFD  BY  ATMOSPHERIC  EXPLOSIONS 

The  coupling  o.'  atmospheric  disturbances  was  considered 
in  the  Seismic  Data  Laboratory  and  the  work  was  carried  over 
into  the  subject  project.  The  computations  were  performed  for 
explosions  over  a  half  space,  both  water  covered  and  solid, 
and  for  explosions  over  an  atoll. 

a.  Rayleigh  waves  excitation  over  a  half  space 

Harkrider  made  available  to  us  his  program  (Harkrider, 
1964a,  1964b)  for  calculating  far-field  modal  propagation  of 
acoustic-gravity  waves  in  an  atmosphere  consisting  of  iso¬ 
thermal  layers  overlying  a  rigid  half  space,  and  a  similar 
program  for  calculating  Rayleigh  wave  modal  propagation  in 
a  layered  halfspace  underlying  a  vacuum.  These  two  programs 
were  merged  in  order  to  model  the  coupling  of  atmospheric 
and  seismic  wave  propagation,  specifically  the  generation  of 
Rayleigh  waves  by  point  sources  in  the  atmosphere.  It  was 
found  that  the  loading  of  the  atmosphere  on  the  solid  earth 
affected  the  Rayleigh  wave  phase  velocity,  group  velocity, 
and  medium  response  by  about  the  ratio  of  atmospheric  sea- 
level  density  to  the  density  of  the  earth's  crust,  i.e., 
about  0.01  percent. 

The  suite  of  programs  used  in  this  project  to  predict 
Rayleigh  wave  excitation  are  described  by  Kerr  (1971).  In  the 
first  program  segment,  the  atmospheric  eigenfunctions  are 
calculated  and  used  to  construct,  the  source-height  excita¬ 
tion  function.  Then  the  actual  far-field  Rayleigh  wave¬ 
forms  were  synthesized  by  the  Aki  algorithm  (Aki,  1960), 
using  the  source  exitation  and  assumptions  about  the  solid 
earth  response  and  Rayleigh  wave  dispersion  characteristics. 
In  the  research  reported  here  we  used  two  basic  earth  models: 
the  Gutenberg  continental  model,  which  has  a  low- velocity 
zone  (Ben-Menahem  and  Harkrider  1964,  p.  2610)  and  the 
Anderson-Toksoz  oceanic  model  (Harkrider  and  Anderson,  1966, 
p.  2970).  Earth  and  atmospheric  layering  characteristics 
were  as  usual  assumed  constant  between  source  and  receiver. 

Calculations  were  made  using  a  wide  variety  of  source 
heights  and  yields.  Several  types  of  sources  were  used: 
mass-injection  sources  (Harkrider  and  Flinn,  1970),  energy- 
injection  sources  (Pierce,  1968),  and  scaled  surface  over¬ 
pressures  (Murphy,  1971).  In  another  report  we  showed 
(Harkrider  et  al.,1972)  that  the  energy- inj ect ion  source 
is  the  most  reasonable  of  those  source  types.  Waveforms 
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and  spectra  were  calculated  for  yields  nf  i  vt  m  ,rr 
burst  heights  of  0.30  to  80  kmT.  Y  1  kT  t0  10  m  at 

These  calculated  data  were  then  ansiv7oH 
source  types,  and  to  find  characteristils^of  the  waveP- ^ 
afndrnyie“dd.SPeCtra  KhlCh  mi*ht  be  ''-gnos tic^o^Vrst6;, eight 
Our  conclusions  are: 

A.  Source  type: 

1.  For  the  generation  of  Ravleieh  wavps 
injection  and  energy- inj ection  sources  are  equivalent." 

_  f  2/.The  Rayleigh  wave  amplitude  for  source  heiah-s  nf 

B.  Earth  model: 

C.  Yield-height  diagnostics: 

yield 

bOTstT'T^  °n  theeyieTdrasfwell  as'on’the"” 

etr|l.h(1972).See  OUr  "10re  detailed  discussion  in  Harkrider 

of  spectral  lEe^ibw^MriSdTt!!  (i-e'Vthe  ratio 
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3  The  spectral  splitting  ratio  spectrum:  This  1S 
simply  *  <R>t  plotted  as  a  function  of  the  splitting  period 
T  These  spectra  are  similar  for  all  yields  at  a  given 
birst  height,  but  have  distinct  shapes  for  different  burst 
heights .  Unfortunately,  the  variation  of  these  sPec”^,J'l‘h 
hnrft  heieht  is  irregular,  and  it  does  not  appear  possible 
M Take  pledictioJis  lor  heights  which  were  not  included  in 
our  calculation  grid. 
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VI.  RAYLEIGH  WAVE  EXCITATION  OVER  AN  ATOLL 


The  Rayleigh  waves  generated  by  an  atmospheric  explosion 
over  an  island  atoll  are  not  easy  to  study  analytically: 
coupling  of  energy  from  the  source  into  the  earth  takes  place 
over  the  shallow  water  in  the  middle  of  the  atoll,  but  a 
transition  takes  place,  well  within  the  critical  distance  for 
Rayleigh  wave  development,  to  the  deep  ocean  whose  properties 
determine  the  spectrum  of  the  waveform  recorded  at  large 
distances.  Approximate  theoretical  coupling  techniques  are 
hence  inapplicable. 

Consequently,  we  undertook  to  study  numerical  approxi¬ 
mations  to  the  equations  which  describe  wave  propagation  in 
inhomogeneous  media.  The  first  step  was  to  use  a  finite 
difference  approximation,  in  which  time  and  space  deriva¬ 
tives  of  the  wave  functions  and  the  material  properties  are 
approximated  by  differences.  After  a  large  amount  of  effort 
and  computer  time,  and  lengthy  discussion  with  mathematicians 
and  seismologists  who  have  been  working  in  this  field  for 
several  years,  it  became  clear  that  the  problem  at  hand  is 
beyond  the  present  state  of  the  art  of  finite  difference 
calculations.  Workers  such  as  Alterman  and  Karal  always 
deal  with  problems  in  which  the  source  time  function  is 
slowly  varying  (with  respect  to  the  time  mesh  size),  and 
the  material  properties  vary  smoothly  and  gradually,  and  it 
turns  out  that  the  reason  is  that  abrupt  material  discon¬ 
tinuities,  and  shock  wave  inputs  --  especially  at 
discontinuous  surfaces  such  as  the  water  surface  --  cannot 
be  handled  satisfactorily  by  thij  method. 

Consequently,  we  turned  to  a  more  sophisticated  pro¬ 
cedure,  the  finite  element  scheme.  In  this  approximation 
the  material  is  divided  into  small  or  large  homogeneous 
polygonal  elements,  which  are  welded  together  along  their 
edges.  The  elastic  wave  equation  is  solved  analytically 
in  each  element,  under  the  single  assumption  that  the 
displacement  components  vary  linearly  along  the  edges  of 
each  element.  The  method  works  because  the  displacements 
at  the  vertices,  as  calculated  within  each  of  the  contiguous 
elements,  must  all  be  the  same.  This  fact  leads  to  condition 
equations  which  can  be  solved  straightforwardly.  A  static 
problem  is  solved  at  each  time  step,  and  the  result  is  used 
as  input  to  the  next  time  step. 

By  a  literature  search  at  DDC,  we  located  a  highly 
sophisticated  finite  element  program  -  SLAM  -  which  had  been 
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developed  for  the  Air  Force  by  the  Illinois  Institute  of 
Technology.  We  obtained  the  program  documentation  --  five 
inches  thick  --  and  studied  the  program  during  a  nine- 
month  delay  waiting  for  the  program  deck  itself.  When  it 
arrived,  we  found  that  it  was  an  older  version  written 
for  the  CDC  6600,  and  contained  many  errors  --  both  logical 
errors  and  elementary  FORTRAN  mistakes.  It  was  clear  that 
what  we  received  was  not  a  working  version  of  the  code, 
and  we  were  unable  to  obtain  a  working  version  despite 
repeated  attempts. 

Part  of  the  program  was  written  in  CDC  6600  assembly 
language,  so  we  obtained  an  in-house  terminal  to  that  type 
of  machine,  and  proceeded  to  correct  the  program  and 
experiment  with  it  on  problems  whose  solutions  are 
analytically  known.  The  total  delay  involved  in  obtaining 
the  program  and  getting  it  working  was  clearly  less  than  would 
have  been  required  to  develop  such  a  program  fiom  scratch. 

Several  modifications  were  made  to  the  program:  we 
made  it  possible  to  vary  the  time-step  acceleration  between 
time  steps;  a  general  implicit  scheme  was  installed,  and 
finally,  a  fundamental  error  involving  the  boundary  at  the 
free  surface  was  detected  and  repaired.  The  last  involved 
the  fact  that  although  the  total  force  on  each  of  the  node 
points  at  the  free  surface  vanished  identically,  the  vertical 
and  radially  tangential  stresses  did  not.  Of  course,  the 
seismic  boundary  condition  is  that  the  stresses  vanish. 

The  results  obtained  within  the  contract  period  were 
disappointing.  One-dimensional  problems  were  solved  correctly, 
but  the  simplest  two-dimensional  problems  (say,  a  line 
source  over  a  homogeneous  halfspace)  gave  answers  which  did 
not  agree  with  analytically  known  solutions.  Arrival  times 
of  P,  S,  and  Rayleigh  waves  were  correct,  but  the  Rayleigh 
wave  was  followed  by  long  ringing  oscillations  which  are 
clearly  artifacts.  The  reason  for  the  ringing  is  not  under¬ 
stood. 

Further  work  on  this  program  and  its  problems  is 
being  continued  under  another  contract  at  the  Pennsylvania 
State  University. 
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